Robot Navigation in Complex Workspaces Employing Harmonic Maps and Adaptive Artificial Potential Fields

In this work, we address the single robot navigation problem within a planar and arbitrarily connected workspace. In particular, we present an algorithm that transforms any static, compact, planar workspace of arbitrary connectedness and shape to a disk, where the navigation problem can be easily solved. Our solution benefits from the fact that it only requires a fine representation of the workspace boundary (i.e., a set of points), which is easily obtained in practice via SLAM. The proposed transformation, combined with a workspace decomposition strategy that reduces the computational complexity, has been exhaustively tested and has shown excellent performance in complex workspaces. A motion control scheme is also provided for the class of non-holonomic robots with unicycle kinematics, which are commonly used in most industrial applications. Moreover, the tuning of the underlying control parameters is rather straightforward as it affects only the shape of the resulted trajectories and not the critical specifications of collision avoidance and convergence to the goal position. Finally, we validate the efficacy of the proposed navigation strategy via extensive simulations and experimental studies.


Introduction
The navigation of autonomous robots in cluttered environments is a widely studied topic in the field of robotics. Popular methodologies that have been employed in the related literature to address it include, but are not limited to, configuration space decomposition approaches [1,2]; probabilistic sampling methods such as rapidly exploring random trees [3,4], probabilistic roadmaps [5,6] and manifold samples [7,8]; and optimal control strategies such as receding horizon control [9,10] and path homotopy invariants [11,12]. Apart from the aforementioned discrete methods regarding the workspace and/or the decision domain, Artificial Potential Fields (APFs) that were originally introduced in [13] generally provide a simpler means of encoding collision avoidance specifications, with their negated gradient functioning as a reference motion direction that drives the robot towards the desired goal configuration. As shown in [14], despite their intuitive nature, this class of controllers suffers unavoidably from the presence of unwanted equilibria induced by the workspace's topology, whose region of attraction may not be trivial. In their seminal work [15], Rimon and Koditschek presented a family of APFs called Navigation Functions (NFs) for point and sphere worlds, as well as a constructive transformation for mapping workspaces cluttered by sequences of star-shaped obstacles into such worlds. However, certain design parameters require tedious tuning to eliminate unwanted local minima and render the transformation a diffeomorphism. In practice, this solution suffers by the fact that the allowable values of the design parameters may cause both the potential and the corresponding transformation to vary too abruptly close to the obstacles (the issue of "disappearing valleys" [15]), thus pushing the trajectories of the robot very close to them. Density functions for remedying such drawbacks or adjustable NFs for relaxing some generally conservative requirements are presented in [16,17]. Additionally, attempts to extend the NF framework directly to non-sphere worlds can be found in [18,19]. Finally, a novel approach based on power diagrams which can be used for designing tune-free vector fields for navigation within convex workspaces is also presented in [20].
Artificial Harmonic Potential Fields (AHPFs) constitute an interesting subclass of APFs, since they are free of unwanted local minima by construction. However, no simple method exists for constructing safe (with respect to obstacle avoidance), harmonic potentials even for simple workspaces. AHPFs suitable for navigation in realistic environments were originally utilized in [21], where computationally expensive numerical techniques were employed to solve the associated Dirichlet and Neumann problems. Several extensions of the aforementioned methodology followed [22,23], addressing issues such as numerical precision and computation, dynamic environments, etc. The panel method was also employed in [24][25][26] to build harmonic potentials to coordinate the motion of single and multiple robots in polygonal environments. In [27,28], well-known closed-form solutions of the uncompressed fluid flow around simple geometries was used in order to safely drive a robot among moving obstacles. Harmonic potential fields have also been used in [29,30] to address the Simultaneous Localization and Mapping problem (SLAM) by coordinating the robot motion in unknown environments. Moreover, a methodology based on the evaluation of the harmonic potential field's streamlines was used in [31,32] for mapping a multiply connected workspace to a disk, collapsing inner obstacles to line segments or arcs. In a recent work [33], the problem of designing closed form harmonic potentials in sphere worlds was addressed by the introduction of a diffeomorphism [34], which allows mapping such workspaces to the euclidean plane with some of its points removed. Finally, extensions of this work addressing topologically complex three-dimensional workspaces or multi-robot scenarios by introducing appropriate constructive workspace transformations can be found in [35,36], respectively.

Contributions
We address the navigation problem for a robot operating within a static, compact, planar workspace of arbitrary connectedness and shape by designing a control law that safely drives the robot to a given goal position from almost any initial feasible configuration. The goal of this work is twofold. (A) To cope with the topological complexity of the workspace, we employed numerical techniques in order to build a transformation that maps the workspace onto a punctured disk and delved into the respective construction in detail. We remark that, although the transformation constructed using this method is an approximation of a harmonic map ideal for navigation, our solution benefits from the fact that it only needs a sufficiently fine polygonal workspace description that can be easily acquired in practice (e.g., through SLAM), contrary to [15,34,36] that require an explicit representation of the workspace boundaries (i.e., as the level sets of sufficiently smooth functions). Moreover, unlike the solutions proposed in [15,36], our approach does not require the decomposition of the workspace obstacles into sequences of simpler overlapping shapes and computes the desired transformation in one step. (B) To steer the robot to its desired configuration, we employed a control law based on closed-form AHPFs coupled with adaptive laws for their parameters to eliminate the necessity of explicitly defined local activation neighborhoods around the workspace boundaries for ensuring collision avoidance. Our approach is reactive (closed loop) since it selects the velocity of the robot based on the positions of the robot, the desired goal and the workspace boundary. As such, it is more robust against position measurement errors than other open loop approaches such as configuration space decomposition approaches [2] or probabilistic sampling methods such as rapidly exploring random trees [4], probabilistic roadmaps [6] and manifold samples [8], where an open loop path is initially extracted and executed by a trajectory tracking controller. In this way, even small position errors risk the safe execution of the calculated plan. We remark that our overall control scheme only requires solving a computationally expensive problem once for a given static workspace, independent of the robot's initial and goal configurations, in contrast to the solutions presented in [21,22]. Finally, we adapt our methodology to the class of differential drive robots, which are commonly encountered in real-world applications and propose an algorithm that decomposes the overall workspace into small neighbouring subsets to render the problem of addressing large workspaces tractable. An overview of the proposed methodology's pros and cons compared to alternative transformations and potential fields can be seen in Tables 1 and 2, respectively. Preliminary results were included in our conference paper [37]. We have to stress though that the algorithmic calculation of the harmonic map is given in the present work, along with a rigorous formulation of the panel method. A modification of the adaptive laws for the parameters of the underlying potential field is also introduced to simplify the tuning process by eliminating the necessity of heuristically defined local activation neighborhoods around the workspace boundaries for ensuring collision avoidance. Moreover, an extension for tackling the navigation problem under unicycle kinematics is also provided. Finally, new comparative simulation results are provided to highlight the strong points of the proposed method with respect to other related works, accompanied by an experiment employing an actual robot navigating within a complex office workspace. Table 1. Comparison between the Harmonic Transformation (HM) proposed in this work and the (i) Star-to-Sphere Transformation (SST) [15], (ii) Multi-Agent Navigation Transformation (MANT) [36] and (iii) the Navigation Transformation (NT) [34]. Although HMs require global knowledge of the workspace's geometry to be constructed, HMs are infinitely differentiable and require the domain to be represented by closed polygonal curves (which can be easily obtained using SLAM methodologies), unlike the alternatives that require the domain boundaries to be represented as sets of sufficiently differentiable implicit equations.

Geometry Representation
Global Analytic  [15], Harmonic Navigation Functions [33] and approximate Harmonic Potential Fields obtained using numerical techniques [21]. Unlike RKNFs that require tuning for ensuring convergence to the goal from almost all initial configurations and HNFs that require tuning for guaranteeing collision avoidance with the workspace boundaries, the proposed control law enjoys both properties by design.

Preliminaries
We use D r (x) to denote an open disk with radius r > 0 centered at x ∈ R 2 . Additionally, D and ∂D denote the closed disk and circle with unit radii centered at the origin of R 2 , respectively. Furthermore, let I N {1, 2, . . . , N} and I N {0} I N . Given sets A, B ⊆ R n , we use cl(A), ∂A, int(A) and A to denote the closure, boundary, interior and complement of A with respect to R n , respectively, and A \ B to denote the complement of B with respect to A. Furthermore, we use 0 N and 1 N to denote the all-zeros and all-ones column vectors of length N, respectively, and 0 N×M to denote the N × M zero matrix. We also define 1 k N×M , k ∈ I M as the N × M matrix whose k-th column is equal to 1 N and every other column is equal to 0 N . Given a vector function f (x), we use ∇ x f to denote its Jacobian matrix. Furthermore, given an arc C, we use |C| to denote its length. We will also say that a set A is attractive (repulsive) under a potential function ψ when there exists a point p 0 / ∈ cl(A) such that if we initialize at p 0 and move along the negated gradient of ψ, we will converge (not converge) to ∂A. Finally, a potential function ψ is called harmonic if it satisfies the Laplace equation, i.e., ∇ 2 ψ = 0, where ∇ 2 denotes the Laplacian operator. An important property of harmonic functions is the principle of superposition, which follows from the linearity of the Laplace equation. Moreover, the extrema of a non-constant harmonic function occur on the boundary of the domain of definition, thus excluding any local minima/maxima within it (a desirable property for motion planning).

Problem Formulation
We consider a robot operating within a compact workspace W ⊂ R 2 bounded by a single outer and a finite set of inner disjoint Jordan curves (a Jordan curve is a non-selfintersecting continuous planar closed curve), which correspond to the boundaries of static obstacles. It is assumed that W can be written as: where W i , i ∈ I N denote regions of R 2 that the robot cannot occupy (see left subplot in Figure 1). Particularly, the complement of W 0 is considered to be a bounded, simply connected region that may also include a strict subset of its own boundary (this corresponds to cases when we wish to place the robot's goal configuration on some part of the workspace outer boundary which is not physically occupied by an actual obstacle, e.g., the door of a compartment (refer to Section 5.2 for more details)) and W 1 , W 2 , . . . , W N are assumed to be closed, simply connected compact sets that are contained in W 0 and are pairwise disjoint. Let p = [x, y] T ∈ R 2 denote the robot's position and assume that the robot's motion is described by the single integrator model:ṗ where u ∈ R 2 is the corresponding control input vector.

Problem 1.
Our goal is to design a control law to successfully drive a robot with kinematics (1) towards a given goal configuration p d ∈ W from almost any feasible initial configuration p init ∈ W, while ensuring collision avoidance, i.e., p(t) ∈ W for all t ≥ 0.

Remark 1.
The results presented in this work can be readily employed for the navigation of disk robots with radius R > 0 by appropriately augmenting the workspace boundaries with the robot's size.

Harmonic Maps for Planar Navigation
In this section, we present a methodology that maps the robot's workspace onto a punctured unit disk, over which the robot's control law is designed. Particularly, our goal is to construct a transformation, T : cl(W) → D, from the closure of the robot's configuration space cl(W) onto the unit disk D with the following properties: 1.
T(·) maps the outer boundary ∂W 0 to the unit circle ∂D; 2.
T(·) maps the boundary ∂W i , i ∈ I N of each obstacle to a distinct point T(·) is a diffeomorphism for all p ∈ int(W).
To that end, we compute a transformationT(p) = [ũ(p),ṽ(p)] T , withũ(p) andṽ(p) being harmonic functions with respect to p, by approximating an ideal harmonic map T that meets the aforementioned properties and the existence of which was proven in Theorem 2 of [38], accompanied by sufficient conditions that render it a diffeomorphism as outlined in the corresponding proof. Particularly, this theorem implies that given an orientation-preserving, weak homeomorphism T ∂ : ∂W 0 → ∂D (such a transformation can be easily obtained for any given planar Jordan curve C by (1) arbitrarily selecting a point p o on C, (2) defining (p), ∀p ∈ C as the length of the arcp o p, assuming one travels from p o to p on C while having the curve's interior to its left and (3) choosing T ∂ (p) = [cos(2π (p)/L), sin(2π (p)/L)] T , where L = |C|) from the workspace outer boundary ∂W 0 to the boundary of the unit disk, the harmonic map T that satisfies the conditions: with n p denoting the unit vector that is normal to the boundary at the point p ∈ ∂W i , i ∈ I N , is a diffeomorphism that maps cl(W) to the target set and collapses each inner obstacle W i onto a distinct point q i within its interior (see Figure 1). Note that the coordinates of q i , i.e., the images of the internal obstacles, are not explicitly specified but are computed as part of the solution, such that the aforementioned constraints are satisfied. Given that closed-form solutions to the aforementioned problem are generally not available for non-trivial domains, in this work, we employed numerical techniques and particularly the Panel Method [24,39,40] (similar formulations can be obtained by employing other numerical techniques such as the Boundary Element Method (BEM), the Finite Element Method (FEM) or the Finite Differences Method (FDM)) in order to construct a harmonic mapT that sufficiently approximates T. As such, by subdividing separately the workspace's outer and inner boundaries intoM 0 ,M 1 , . . . ,M N number of elements (see Figure 2), we define the components ofT(p) = [ũ(p),ṽ(p)] T as follows: whereM C is the number of control parameters per element,Ẽ ij denotes the j-th element of the i-th boundary's approximation,p i,j (s) : [0, |Ẽ ij |] →Ẽ ij is a bijective parameterization of E ij ,G ijl : [0, |Ẽ ij |] → R is the shape function corresponding to the l-th control parameter ofẼ ij andC x ijl ,C y ijl ∈ R are control parameters that need to be appropriately selected so thatT satisfies properties 1-3 for all l ∈ IM C , j ∈ IM i and i ∈ I N . It is worth noting that for common choices ofG ijl (e.g., constant or linear shape functions) and simple types of E ij (e.g., line segments), the integral in (4) can be easily evaluated to obtain a closed-form expression forH ijl . As an illustration, for a line segment elementẼ ij with two control parameters (i.e.,M C = 2), a typical choice for linear shape functions (see Figure 2) is s/|Ẽ ij | for the corresponding parameterization, wherep i,j,A ,p i,j,B are the element's end-points. To obtain the unknown control parameters as well as the images of the workspace's inner obstacles, one needs to solve two independent linear systems of equations: for the unknown vectors: The matrixÃ and the right hand side vectorsB x andB y are constructed by selecting a set of ∑ i∈I Nm i arbitrary pointsp i,j (a typical strategy is to select the pointsp i,j uniformly on the outer and inner boundaries of the given domain) such that a)p i,j ∈ ∂W i for all j ∈ Im i and i ∈ I N and b) ∑ i∈I Nm i =M C ∑ i∈I NM i , and evaluating (2) and (3) at those points as follows: Notice that by discretizing the workspace boundaries into a large number of sufficiently small elements, the overall approximation error between the solutionT of the aforementioned linear problem and the exact transformation T can be rendered arbitrarily small (see [39,40]). However, the complexity of constructing the mapping is of order O(M 3 ), whereM denotes the number of total elements of the mapping (i.e., the complexity of the solution of the dense system of linear Equation (5)). Nevertheless, the construction of the transformation, which is the main computational bottleneck, is performed only once at the beginning. Additionally, apart from the straightforward user-defined homeomorphism T ∂ on the workspace boundary, no tedious trial and error tuning is needed to extract the diffeomorphic transformationT, in contrast to other related works such as the Star-to-Sphere Transformation (SST) [15], the Multi-Agent Navigation Transformation (MANT) [36] and the Navigation Transformation (NT) [34]. Discretization of a given domain's boundary using line segment elements. By convention, the outer boundary is considered to be clockwise oriented, whereas inner boundaries are counterclockwise oriented. The normal direction of each element is depicted using green colored vectors. Furthermore, the values of the two linear shape functionsG 0,3,1 andG 0,3,2 are plotted along the associated elementẼ 0,3 .

Control Design
To address Problem 1, we equip the robot with the aforementioned transformation q = T(p) from the closure of its configuration space W onto the unit disk D and an artificial potential ψ(q, k) augmented with an adaptive control lawk = f k (q, k) for its parameters k = [k d , k 1 , k 2 , . . . , k N ] T . The robot velocity control law is calculated as follows: where J(p) denotes the Jacobian matrix of T(p), s(q, k) ≥ 0 is a continuously differentiable gain function given by: with and K u , γ, p and v are scalar constants such that K u , v > 0 and γ, p ∈ (0, 1). More specifically, s(q, k) consists of two individual terms, with the first vanishing as the robot approaches the workspace's outer boundary (and its distance from the unit circle is less than p ) and the second vanishing when the robot's velocity points away from the disk's center. The scalar parameter γ can be used for adjusting the contribution of each respective term of s(q, k). Finally, ψ is a harmonic artificial potential field defined on the image T(W) of the workspace W and whose negated gradient −∇ q ψ(q, k) defines the direction of the robot's motion in the real workspace W via the inverse Jacobian J −1 (p). By design, the resultant vector field precludes collisions between the robot and the workspace's inner obstacles and renders the goal configuration almost globally attractive except for a set of measure zero initial configurations. However, since W 0 may not be repulsive under ψ for an arbitrary, fixed selection of k, we also introduce the adaptive law f k (q, k) which, along with s(q, k), guarantees forward invariance of the workspace without compromising the convergence and stability properties of the overall system. The following subsections elaborate on each component of the proposed control law individually.

Artificial Harmonic Potential Fields
We construct an artificial harmonic potential field on the disk space D employing point sources placed at the desired configuration q d = T(p d ) as well as at the points q i = T(∂W i ), ∀i ∈ I N that correspond to the inner obstacles, as follows: where k d > 0 and k i ≥ 0 denote harmonic source strengths which vary according to adaptive laws that are presented later. An interesting property of the above potential field, which stems from the maximum principle for harmonic functions, is that, for fixed k, the only minima of φ are located at q d and, possibly, at infinity. As a direct consequence of this property, the Hessian ∇ 2 q φ computed at a non-degenerate critical point of φ in our domain's interior has one positive and one negative eigenvalue with the same magnitude, e.g., λ and −λ with λ > 0.
Next, we define a reference potential ψ based on φ, which is given by: where w φ is a positive scaling constant. Note that ψ maps the extended real line to the closed interval [0, 1]. As tanh φ/w φ is a strictly increasing function, the only critical points of ψ are the ones inherited from φ with their indices preserved. Furthermore, the gradient of ψ with respect to q, given by is well defined and bounded for all q ∈ D.
If the workspace was radially unbounded, selecting k fixed with k d > ∑ N i=1 k i would render the potential field (10) sufficient for navigation. The author of [33] addresses bounded workspaces that are diffeomorphic to sphere worlds by simply mapping the outer bounding circle to infinity. In this work, we would like to be able to place q d on regions of ∂D that are not physically occupied by obstacles (such as passages to other compartments, see, for example, Section 5.2); thus, we cannot follow the same procedure since that would render the effect of the sole attractor on the robot null. Instead, we design appropriate adaptive laws for the parameters k of φ to render the outer boundary repulsive and establish the forward completeness of the proposed scheme at all times.
Before proceeding with the definition of the adaptive law, we first state two propositions that will be used in the subsequent analysis, the proofs of which can be found in the Appendix A.
Proposition 2. If k i are non-negative and bounded, there exists k d > 0 such that ψ is Morse for all k d ≥ k d .

Adaptive Laws
We now present the adaptive lawk = f k (q, k) that updates the parameters of the potential field ψ. Its primary goal is to render (a) the workspace outer boundary repulsive and (b) any critical point of φ in the vicinity of the robot non-degenerate, a property that will be used later in the analysis. In particular, we consider f k of the form: where w i and g i , i ∈ I N , as well as h i , i ∈ I N , are functions to be defined later, k i , i ∈ I N are desired upper bounds for k i , λ denotes the non-negative eigenvalue of ∇ 2 q φ, K k is a positive control gain and 1 and 2 are small positive constants. The continuously differentiable switch ξ 1 (x; ) and functions i (q) are, respectively, given by: According to Proposition 1, our first requirement can be accomplished by designing f k to reduce k i as the robot approaches ∂D. To do so without compromising the inherent inner obstacle collision avoidance properties of φ, we need to also ensure that each k i does not vanish within some neighborhood of q i for all i ∈ I N . To that end, firstly we define g i , employing the smoothly vanishing function defined in (8) to serve as pseudo-metrics of the alignment between the robot's velocity and the directions towards the goal and inner obstacles, respectively, given by: with where α ∈ (0, 1] is a fixed constant that is used for selecting the desired alignment between the robot's motion and the direction to the goal. We also define the accompanying weights w i as follows to ensure that only one term of (12) dominates as the robot approaches a particular boundary of W: with for a scalar constant 3 ∈ (0, 1) in (15) and some integer m < −1 in (16) that serves as a smooth under-approximation of min j =i (r j ), i ∈ I N . Finally, the weights h i , i ∈ I N are defined as follows: whose purpose is to accelerate the decay of those k i that contribute the most to the component of ∇ q ψ that pushes the robot toward the workspace's outer boundary.
Regarding the second requirement, as shown in Proposition 1, selecting a k d above a certain threshold is sufficient to render φ free of degenerate equilibria. On the other hand, for a given k i , increasing k d steers the robot closer to the workspace's inner obstacles. Nevertheless, since the robot may never actually enter the vicinity of a degenerate equilibrium, instead of setting k d sufficiently large a priori, the adaptive law for the parameter k d is introduced to increase k d only when it is actually needed, thus alleviating the aforementioned shortcoming.

Stability Analysis
Let us consider the overall system:ż where z = (q, k) and , Ω = T(W). Note that Ω consists of int(D), possibly with a subset of ∂D, with the points q i removed. In this section, we elaborate on the stability properties of (17) under the proposed control scheme (6) and (12). First, we formalize the safety properties of the closed-loop system dynamics, which guarantee that our robot does not collide with any obstacle.
Proof. For the proof, refer to the Appendix A.
Having eliminated the possibility of the robot colliding with the workspace's boundaries, we proceed by showing that all critical points of ψ, where (17) may converge to, are either non-degenerate saddles or q d . Additionally, we show that the latter is a stable equilibrium. Proposition 4. The artificial potential ψ decreases along the trajectories of the closed-loop system and its time derivative vanishes only at its critical points. Additionally, the preimage of q d is a set of stable equilibria of (1).
Proof. For the proof, refer to the Appendix A.
Proposition 5. Let z = (q , k ) be a critical point of the closed-loop system dynamics with q ∈ Ω \ {q d }. Then, q is a non-degenerate saddle point of ψ.
Proof. For the proof, refer to the Appendix A.
Finally, we conclude this section with the main theoretical findings.  (6) and (12) converges safely to p d , for almost all initial configurations, thus addressing successfully Problem 1.
Proof. For the proof, refer to the Appendix A.

Remark 2.
Owing to the adaptive laws (12) that modify the harmonic source strengths online to secure the safety and convergence properties at all times, the selection of the fixed control parameters in the proposed scheme, i.e., K u , γ, p , v , w φ , K k , 1 , 2 , α and 3 , is straightforward as it affects only the trajectory evolution within the workspace and not the aforementioned critical properties. Consequently, their values should be set freely as opposed to NFs, where the selection of the main parameters severely affects the convergence properties of the adopted scheme and cannot be conducted constructively for generic workspaces of arbitrary topology.

Extensions
In this section, we present certain extensions of the proposed approach to (a) address the safe navigation problem for unicycle robots which are frequently encountered in many application domains and (b) tackle computational complexity issues that affect the numerical computation of the harmonic map presented in Section 3 as the size of the workspace increases.

Unicycle Robot Kinematics
In this subsection, we consider robots whose motion is subjected to Pfaffian constraints of the form: where θ ∈ [0, 2π) denotes the robot's orientation, n θ = [cos(θ), sin(θ)] T , and v, ω ∈ R are control inputs corresponding to the robot's linear and angular velocities, respectively. First, let us define the robot's kinematics in the image of the configuration space via the proposed transformation as follows:q = nθvθ =ω .
Note that the orientations θ andθ are related by: To safely drive the robot to its goal configuration, we consider the following control laws:v with K v , K ω ∈ R positive constant gains, n ⊥ θ = [− sin(θ), cos(θ)] T and Additionally, we need to employ a modified version of the adaptive law for the potential field parameters, which is obtained by substituting s with s v in (12) and (13) and respectively, in (14). Finally, by expressing the aforementioned control laws to the robot's actual configuration space, we obtain: where ω dq and ω dθ are terms corresponding to angular velocities induced by translational and rotational motion of the robot in the workspace's image, respectively, given by: J p denoting the directional derivative of J p along n θ . The stability properties of the aforementioned closed-loop system dynamics are formalized below.

Theorem 2.
The workspace W is invariant under the dynamics of (18) equipped with the proposed control law. Additionally, the robot will asymptotically converge either to an interior critical point of φ or to the pre-image of q d , which is stable.
Proof. For the proof, refer to the Appendix A.
Remark 3. The result of Theorem 2 is weaker compared to that of Theorem 1, since there is no guarantee that the set of configurations which converge to a critical point of φ (other than the pre-image of q d ) has Lebesgue measure zero.

Atlas of Harmonic Maps
As the size of the workspace increases, the problem of computing the transformation T grows in complexity as well, because the resources required by commonly employed numerical techniques that can solve the problem presented in Section 3 are polynomial in the number of elements used for representing W. Alternatively, to cope with large workspaces efficiently, we propose instead the construction of an atlas A {(P i , T i ) | i ∈ I N A } obtained by separating the workspace W into N A overlapping subsets P i ⊂ W, such that i∈I N A P i = W and constructing a separate harmonic map T i for each P i (see Figure 3). This essentially allows us to solve many small (and computationally less intensive) problems instead of a large one, thus reducing the overall resources required for addressing a given workspace. Therefore, given such a partitioning of W, we define the graph G = (V, E ), where V = {P i | i ∈ I N A } denotes the set of corresponding nodes (workspace partitions) and E denotes the set of edges between the elements of V, with each edge indicating a feasible transition from one partition to another, i.e., (i, j) ∈ E if and only if Ä cl(P i ) ∩ cl Ä P j ää = ∅. Note that G is undirected by definition, i.e., (i, j) ∈ E only if (j, i) ∈ E . Additionally, since the workspace is connected, G should also be connected. Thus, for a given atlas A, an initial configuration p init and a final configuration p d , we can employ standard graph search algorithms to obtain a sequence of indices S = {s 1 , s 2 , . . . s n } corresponding to partitions that the robot can tranverse to reach its goal. (In general, more than one such sequence of partitions may exist connecting the initial and the final configurations. However, the selection of one that corresponds to some sort of "optimal" path is beyond the scope of this work.) Additionally, note that since the partitioning of W does not need to be fine, the size of G will generally be small, rendering the cost of finding S negligible.
We now concentrate on how the transition between two consecutive elements of S is implemented. Let C i,j cl(P i ) ∩ cl Ä P j ä denote the common region of cl(P i ) and cl Ä P j ä and let B i,j ∂P i ∩ P j denote the set of points on the boundary of P i that also belong to P j and are not occupied by obstacles for all i ∈ I N A and all j such that (i, j) ∈ E . Without loss of generality, we assume that A is constructed such that the sets B ,i ∩ B ,j are either empty or consist of isolated points. We note that in order to successfully complete the transition between two consecutive nodes P i and P j of S, it suffices for the robot to reach any single point of B i,j from P i . We also observe that each B i,j may consist of one or more disjoint components B i,j , ∈ L(i, j), with L(i, j) being some valid indexing of those. By exploiting the fact that Theorem 2 [38] imposes a weak homeomorphism requirement on T i , we can construct each T i such that each disjoint subset of ∂P i collapses into a separate point, i.e., T i (B i,j ) = q i,j ∈ ∂D (see Figure 3), which, in turn, implies that selecting q i,j as an intermediate goal configuration suffices to render the entire B i,j attractive. Building upon this fact, for each consecutive pair of P i and P j in S, we (arbitrarily) select a B i,j and we construct a transformation T i : P i → D, with q [i] = T i (p), and artificial potential field φ i (q [i] , k [i] ) with goal configuration q [i] d = q i,j . Additionally, to smooth the transition between consecutive partitions, when they overlap, we propose the following modified control law for the robot: where u [i] and u [j] denote the control inputs as defined in (6) and evaluated using ψ i , T i and ψ j , T j , respectively; the function η t,i,j : C i,j → [0, 1] is any smooth bump function such that ; and 4 > 0 is a fixed parameter. What this modification essentially does is incrementally add an extra component, with the direction of ∇ p ψ j , to the robot's velocity when that component is cosine similar (two vectors u and v are cosine similar if their inner product is positive) with ∇ p ψ i . We note that η c,i,j → 1 and η t,i,j → 1 as the robot approaches the boundary of the corresponding partition. We also remark that once the robot has completed its transition to P j , we do not concern ourselves with u [i] anymore, i.e., u = u [j] even if p returns to C i,j . The overall scheme employed for navigating a holonomic robot to its goal configuration using an altas constructed as described above can be found in Algorithm 1.
Regarding the stability analysis of the modified system, by following the same procedure as in Section 4.3 and by virtue of η c,i,j , it is trivial to verify the following statement. Theorem 3. System (1) equipped with Algorithm 1 converges safely to a given goal configuration p d ∈ W from almost all initial configurations p init ∈ W.
Proof. For the proof, refer to the Appendix A.

Algorithm 1 Altas-based motion planning scheme for a holonomic robot
Require: A, p init , p d S ←FINDPATHTOGOAL(G, p init , p d ) Initialize k [s] for all s ∈ S. for all i in I n−1 do s, s ← s i , s i+1 Select (arbitrary) such that ∈ L(s, s ). Place goal configuration of ψ s at q s,s . end for Place goal configuration of ψ s n at T s n (p d ). ← 1 loop if = n or p ∈ P s \ P s +1 then Update p using (6) and k [s ] using (12). else if p ∈ C s ,s +1 then Update p using (22)

Simulations and Experimental Results
In order to demonstrate the efficacy of the proposed control scheme, we have conducted various simulation and experimental studies, the results of which are presented in this section. The algorithm that computes the harmonic transformation and its Jacobian was implemented in C++, while the proposed control protocols were implemented in Python. Code implementations can be accessed at https://github.com/maxchaos/hntf2d (accessed on 16 April 2023). All simulations were carried out on a PC with an Intel i5 processor operating at 2.2 Ghz, with 4 GB RAM and running a GNU/Linux operating system. For more details regarding both simulations and experiments, the reader may refer to the accompanying video material at https://youtu.be/I6WUS81iDh4 (accessed on 16 April 2023).

Simulations-Full Workspace Transformation
In the first case study, a single transformation of the entire 8 m × 5 m workspace (see Figure 3) was constructed and the robot was instructed to navigate to various goal configurations starting from the same initial position. The initial configuration and the parameters of our controller were selected such as to better demonstrate the guaranteed collision avoidance properties of our scheme. Particularly, the initial values for the parameters of the adaptive law were selected as k d = 20, k i = 1 and k i = 20 for all i ∈ I 10 . The values of the remaining parameters were K u = 100, w φ = 20, K k = 100, α = 1, p = 0.025, v = 0.1, γ = 0.7, 1 = 0.01, 2 = 0.1 and 3 = 0.1. The goal configurations and the trajectories executed by the robot, both in the real and transformed workspace, are illustrated in Figure 4.
The simulations were conducted using the Euler method with 10 ms steps. Regarding the computational complexity of the control scheme, the construction of the harmonic transformation for this large workspace that was carried out offline once required 5.4 s to complete for a sufficient approximation of the workspace boundary with 3680 segments. Finally, the online computation of the transformation T(p) and its Jacobian J(p) required an average of 6.0 ms per step.

Simulations-Atlas of Harmonic Maps
In this case study, we decomposed the aforementioned workspace into separate partitions (see Figure 3) and constructed a harmonic transformation T i for each one (we adopted the door of each room as the common boundary between neighboring partitions). The robot was initialized at the same position as the previous study and it was instructed to navigate towards the same set of individual goal configurations. The initial values selected for the parameters of the adaptive law were denotes the amount of obstacles inside the corresponding partition. All remaining control parameters were selected as in Section 6.1. The trajectories of the robot are depicted in Figure 5. The time spent to construct the corresponding harmonic transformations varied from 0.019s to 0.211s (depending on the amount of elements required for sufficiently approximating each room, ranging between 320 and 1000 segments) and was significantly much less than the full map construction of the previous case (5.4 s). Additionally, the online computation of T i (p) and J i (p) in each of these rooms required an average time between 1.0 ms and 2.2 ms per step, respectively. Finally, it should be noted that in this case, the workspace inner obstacles were mapped to points further away from the boundaries of the partitions, which is an interesting result as it alleviates possible numerical issues that may arise in the computation of the transformation near the obstacles (the condition number of the Jacobian of the transformation is improved). It should be stressed that the length of the paths in the second case was less (improvement of 0.5 m on average), owing to the fact that the robot gets closer to the workspace boundary since the individual transformations in each room obtain a better conditioned Jacobian (condition number 0.212 against 0.093) and thus are more fine than the first approach, where a transformation is built for the whole workspace.

Comparative Study-Workspace Transformation
In this subsection, we provide a comparative study of the harmonic map presented in this work against readily available workspace transformation methods employed in the motion planning literature. Particularly, we consider four 4 m × 4 m compact workspaces, each associated with a pair of initial and goal positions, and construct appropriate transformations for each one by employing the methodology presented in this work (HM), as well as (i) the Star-to-Sphere Transformation (SST) [15], (ii) the Multi-Agent Navigation Transformation (MANT) [36] and (iii) the Navigation Transformation (NT) [34] (with the aforementioned Star-to-Sphere transformation serving as the underlying map). The trajectories of the robot executed while tracing the line segment connecting the initial and goal configurations in the images of each domain can be seen in Figure 6. We note that manual tuning of the compared transformations was necessary in order to render each a diffeomorphism but without making them too steep around the obstacles. Furthermore, the domain boundaries considered here had to be sufficiently smooth in order for methodologies such as MANT to be applicable. Finally, we remark that the trajectories corresponding to the proposed transformation are, in general, less abrupt compared to the rest, a property attributed to the fact that our approach is global as opposed to the other transformations, i.e., the distortion caused by each obstacle is not limited to some narrow neighborhood around it. The total length, maximum curvature and distance from the obstacles of each executed trajectory can be seen in Tables 3-5, respectively. We can see from these values that the actual trajectories yielded using harmonic maps are among the shorter and smoother ones, although they tend to approach the obstacles more than the rest.

Comparative Study-Control Law
In this subsection, we provide a comparative study of our control scheme against other motion planning methodologies.

APF-Based Schemes
To demonstrate the efficacy of the proposed control scheme compared to other APFbased schemes, we considered the 12 m × 16 m workspace depicted in Figure 7, for which we constructed a harmonic map as described in Section 3. Next, we equipped a holonomic robot with three alternative control laws and instructed it to visit four distinct goal positions using these controllers, starting each time from a fixed initial configuration. Particularly, we considered a conventional navigation function-based controller (NF) [15] augmented by [17], for the selection of its notorious parameter, and a harmonic navigation function-based controller (HNF) [33], in addition to our adaptive control scheme (AHNF) described in Section 4. We note that all three control laws considered here make use of the same underlying harmonic map T constructed as described above in order to drive the robot to its instructed goal positions. The trajectories executed by the robot can be seen in Figure 7. We remark that, in general, our approach steers the robot away from inner obstacles that lie between its initial and goal configurations, unlike "greedy" schemes such as the conventional NF-based controller, while keeping the traced paths shorter compared to HNFs with fixed source weights, a property attributed to the proposed adaptive laws (12) which penalize misalignment between the robot's velocity and the direction towards the goal configuration.  Table 3. Trajectory lengths (m) executed by employing the four alternative transformations in each workspace displayed in Figure 6.  Table 4. Maximum value of curvature (m −1 ) associated with each trajectory displayed in Figure 6.  Table 5. Minimum distance (m) between each robot trajectory and the corresponding workspace boundaries displayed on Figure 6. The total length and distance from the obstacles of each executed trajectory can be seen in Tables 6 and 7, respectively. First, we have to stress that the length trajectory corresponds to the travelled path towards the goal configuration and thus needs to be small, whereas the minimum distance to the workspace boundary refers to the closest point of the trajectory to the workspace boundary and thus needs to be large to have a safe trajectory. Consequently, note from Table 6 that the NF scheme yielded shorter path lengths than the proposed method in two cases (blue and yellow); nevertheless, such paths approach closer to the workspace boundary as indicated in Table 7, thus resulting in more risky paths. On the other hand, the Adaptive Harmonic Potential Field yields a good trade-off between path length and minimum distance to the boundary, since it achieves the shortest paths for two cases without compromising them, as is the case with the NF. On the other hand, the HPF tend to travel around the obstacle closer to the outer workspace boundary and hence exhibit more safe trajectories but they are significantly longer than the other two schemes.

Sampling-Based Scheme
To compare the control scheme proposed in this work against sampling-based methods, we considered a holonomic point-sized robot positioned inside a 6 m × 8 m compact workspace and a desired goal configuration. To complete this task, we employed two different controllers, namely the one proposed in this work and an admissible planner based on an improved probabilistic roadmap method (PRM) [6]. The trajectories executed by the robot using our control law as well as two of the trajectories generated by the PRM-based planner can be seen in Figure 8. The construction of the associated transformation took 31 s to complete for a given boundary approximation made of 7842 elements, whereas the PRM-based planner required approximately 24 s on average over 10 successful runs to yield a solution (we have to stress that we ran 14 trials to get 10 solutions, since four runs did not complete as they exceeded the 500 s calculation time), using the same boundary approximation for collision checking. The robot trajectories exhibited similar lengths in both algorithms (22.5 m for our method against 21.8 m on average), although no path optimization was employed in our case. Additionally, the proposed scheme resulted in a smoother robot trajectory (based on the resulting sequence of points in both cases, we calculated the minimum curvature radius as 0.23 m for our method against 0.12 m on average for the PRM method). On the other hand, note that our approach solves the motion planning problem for any pair of initial and final configurations within the workspace, whereas the sampling-based scheme considers only one go-to problem. Thus, a different initial or final configuration would require a new solution with the PRM method. On the contrary, the proposed transformation needs to be calculated only once to solve the motion planning problem for any pair of initial and final configurations. Finally, it should be noted that for a narrower corridor in Figure 8, the sampling-based approach failed to derive a solution with a reasonable execution time (no solution was calculated within 500 s), since the probability of sampling connected points within this snaky passage reduces drastically. On the contrary, the proposed transformation took 38 s to complete for the same number of elements (i.e., 7842 elements). Figure 7. Trajectories of the robot navigating to four distinct goal configurations (black crosses) with red, green, yellow and blue color starting from the same initial position (black circle) while using various alternative APF-based controllers.  Table 7. Minimum distance (m) between the corresponding workspace boundaries and each trajectory displayed in Figure 7.  8. Trajectories of the robot navigating to its goal configuration (black cross) generated using the proposed control law and a PRM-based planner.

Experiments
In order to verify the results presented in Section 5.1, real experiments were conducted on a non-holonomic robotic platform (Robotnik Summit-XL) operating within the 10 m × 25 m compact workspace that is depicted in Figures 9 and 10. The boundaries of the workspace were obtained using readily available SLAM algorithms and were later augmented with the robot's shape (approximated by a disk). The workspace was partitioned into six overlapping subsets and the robot was instructed to visit three different goal configurations, each located in a different room. An off-the-shelf localization algorithm was employed for estimating the robot's position and orientation using its on-board sensors (laser scanners and RBG-D cameras), providing feedback at approximately 5 Hz to the robot's linear and angular commanded velocities. The construction of the associated transformations over the six subsets of the workspace took from 1.3 s for the simple and smaller partitions with 800 elements to 3.1 s for the more complex ones employing 1500 elements. On the other hand, the evaluation of the mapping as well as its Jacobian took less than 6 ms on average, which was satisfactory given the low position update rate. Note that our algorithm successfully managed to drive the robot safely (the minimum distance to the workspace boundary was 0.15 m when passing through the doors) to its specified goal configurations, as one can verify from the trajectories (see Figure 9, Figure 11 and the accompanying video material). However, an issue that needs to be pointed out is the oscillating behavior that the robot exhibited in the configuration space's image, particularly in subsets p1 and p2 as depicted in Figure 11. Such behavior is attributed both to (a) the relative slow update of the robot's pose estimation and (b) the inversion of the Jacobian which is ill-conditioned close to extremely narrow passages of the domain. Nevertheless, such shortcomings can be alleviated by a better choice of partitions, e.g., by partitioning the domain into more subsets with less complex shapes. As a future research direction, we shall investigate whether the condition number of the Jacobian of the transformation is a fine criterion, since the condition number is usually used to measure how sensitive a function is to changes or errors in the input, and the output error results from an error in the input via the Jacobian.

Conclusions and Future Work
In this work, we employed harmonic map theory to devise a transformation of complex workspaces directly to point worlds that are appropriate for robot navigation. Subsequently, we presented a novel motion planning control scheme based on closed-form harmonic potential fields equipped with appropriate adaptive laws for their parameters, which can safely navigate a robot to its goal state from almost all initial configurations. Additionally, we extended our approach to accommodate the navigation problem of non-holonomic robots and kept the numeric computations tractable for large workspaces.
Regarding future directions, our aim is first to increase the applicability of the proposed navigation framework by addressing partially known dynamic workspaces, which is far from being straightforward. To remedy the issue of calculation time in this case, we shall adopt a sensitivity analysis approach so that we do not solve the whole problem from scratch, but find how the solution deviates when a small change in the workspace occurs. In this way, we envision a reasonable calculation time (except from the first calculation) that would result in an almost real-time calculation of the transformation and thus allow us to consider even moving obstacles in dynamic environments. However, critical issues have to be studied concerning cases where the workspace changes topologically (e.g., in the case of antagonistically moving obstacles) and this results in significant changes in the transformation. In the same vein, switching in the transformation output might raise practical issues such as chattering that have to be carefully considered. Note that the aforementioned research direction could also serve as a first step towards the solution of the multi-robot motion planning problem, where for each robot all other robots should be considered as moving obstacles, operating antagonistically to achieve their goal configurations. Finally, another challenging research direction concerns the extension to 3D workspaces. Unfortunately, the harmonic maps have been studied only for 2D workspaces, since they rely heavily on complex analyses. Nevertheless, we propose to decompose the 3D motion planning problem into several 2D sub-problems, where the proposed solution works, and then combine them (e.g., decompose the motion along the z-axis and on the x-y plane).
Author Contributions: Methodology, P.V. and C.P.B.; Validation, P.V.; Formal analysis, P.V. and C.P.B.; Writing-original draft, P.V.;Writing-review & editing, C.P.B. and K.J.K.; Supervision, C.P.B. and K.J.K. All authors have read and agreed to the published version of the manuscript.
The gradient of φ with respect to q is given by Computing the inner product of ∇ q φ and q yields: Given that all q i lie within int(D), the second term on the right-hand side of (A2) is finite for all q ∈ ∂D. Similarly, the first term on the right-hand side of (A2) is positive for all q = q d . Let q ∈ ∂D \ {q d }. Additionally, the continuity of 1 − q T d q / q − q d 2 and (1 − tanh φ/w φ 2 )/(2w φ ) implies that there exists a closed neighborhood F (q ) of q , not containing q d , where both are positive. Hence, selecting 1 q−q i é ensures that (∇ q φ) T q > 0 for all q ∈ F (q ). Moreover, computing the derivative of d = 1 − q 2 with respect to time for all q ∈ F (q ) and assuming k i < k , ∀i ∈ I N yieldṡ d = 2K u s∇ q ψ T q > 0; thus, the distance from the workspace boundary increases, which concludes the proof.

Appendix A.2. Proof of Proposition 2
Similarly to the proof of Proposition 3 in [33], we proceed by definingq d q − q d , q i q − q i for all i ∈ I N . Let alsoq d q d / q d andq i q i / q i . Accordingly, the Hessian of φ can be computed by: Note that at a critical point of φ it holds that: Substituting (A4) into (A3) and re-arranging the terms yields: Next, we argue that for any given set of radii ρ i > 0 such that D ρ i (q i ), i ∈ I N are disjoint disks that lie entirely within our domain, there exists k d > 0 such that no critical point of φ exists within D \ i∈I N D ρ i (q i ) for all k d > k d . This implies that, by choosing a sufficiently large k d , each critical point of φ belongs to a single D ρ i (q i ). Let q be a critical point and = argmin i∈I N q − q i . To show that ∇ 2 q φ(q ) is not degenerate, it suffices to show that its eigenvalue λ(q ) is positive. We recall that λ is lower bounded by the quadratic formx T ∇ 2 q φx for all x = 1. By considering the direction ofq and after some tedious calculations, we obtain: The first right-hand side term of (A5) is strictly positive. Since all k i are bounded and non-negative, choosing a sufficiently large k d renders the second and third right-hand side terms non-negative. Furthermore, note that the fourth and fifth right-hand side terms are bounded for all q ∈ D ρ (q ). Thus, by choosing a sufficiently large k d , the first three terms of (A5) can be made dominant, thus renderingq T ∇ 2 q φq positive at q , which concludes the proof. Firstly, we will show that the robot cannot escape through the workspace's outer boundary. Let us assume that q → q ∈ ∂D \ {q d }. Then,q → 0 by virtue of (7), since s(q, k) = 0 for all q = 1 with ∇ q φ T q ≤ 0. Additionally, w 0 → 1 and w i → 0, for all i ∈ I N . Thus,k i < 0 holds within a neighborhood of ∂D, while k i > 0, which implies that k i → 0 for all i ∈ I N . Moreover, Proposition 1 dictates that there exists k > 0 for which any point in ∂D \ {q d } is repulsive under ψ. Since (12) dictates that all k i become less than k in finite time, this contradicts our supposition.
Next, we consider collision avoidance between the robot and the inner obstacles. Let us assume that the robot approaches obstacle i. By construction, w i → 1 while ∇ q ψ → 0 and w j → 0 for all j ∈ I N \ {i}. Note that there exists a neighborhood N i of q i such that w 0 = 0 for all q ∈ N i due to continuity of w 0 and ξ 2 (w 0 ; 3 ). Additionally, since the robot is assumed to approach q i ,q T (q − q i ) cannot be identically zero inside N i . As such, as long as k i < k i ,k i ≥ 0 inside N i withoutk i = 0 for all q ∈ N i . This implies that k i → 0 as q → q i , thus rendering q i a local maximum of ψ. Thus, there exists a neighborhood of q i inside which (∇ q ψ) T (q − q i ) > 0, which contradicts our assumption. since ∇ q ψ(q ) = 0. Furthermore, by construction of the adaptive law (12), the Jacobian of f k with respect to z at z is 0 (1+N)×(3+N) . Thus, linearization of the system f z at z yields ô .
Since the top-left block ∇ 2 q φ is invertible at z , using the well-known property of block matrix determinants, we can see that ∇ z f z has two non-zero eigenvalues, particularly the eigenvalues of ∇ 2 q ψ and a zero eigenvalue with multiplicity 1 + N. Thus, ∇ z f z (z ) has exactly one positive eigenvalue, rendering z a saddle point of (17) (Theorem 3.7 [41]).
Appendix A.6. Proof of Theorem 1 In Proposition 4, we have proven thatψ < 0 for all q ∈ Ω \ {q d }, except for the critical points of φ that lie in it. Lasalle's Invariance Theorem (Theorem 3.4 [41]) dictates that system (17) will converge to either (a) the desired configuration q d , (b) the obstacles q i or (c) a critical point z = (q , k ) with q ∈ Ω \ {q d }. We know from Proposition 3 that the critical points of case (b) are repulsive; therefore, no trajectory of the system may converge to them. Regarding the critical point z corresponding to case (c), Proposition 5 dictates that it must be a non-isolated, degenerate equilibrium of the whole of system (17), since ∇ z f z has one positive, one negative and several zero eigenvalues. Let k d be the upper bound of k d that the closed-loop system can possibly attain, as indicated by Proposition 2. In order to prove that the set of initial conditions leading to these points has zero Lebesgue measure, we will study the properties of the gradient-like system (by definition, a gradientlike system is a pair of a scalar cost functions and a dynamical system for which each non-equilibrium initial condition moves the state towards a new one whose cost is less than that of the initial state) (ψ(z), F z,τ (z)) in the domain S z , where the scalar potential ψ(z) is treated as a function to be minimized, the map F z,τ (z) : S z → R N+3 is given by F z,τ (z(t)) z(t + τ) = z(t) + t+τ t f z (z(s))ds for any τ > 0 and S z D × [1, k d ] × [0, k 1 ] × . . . [0, k N ]. Note that S z is convex and closed. Additionally, the map F z,τ (z) is a locally Lipschitz diffeomorphism in S z and S z is forward invariant under F z,τ (z) (by virtue of Proposition 3 and design of adaptive law (12)) for all τ > 0. Furthermore, the unwanted equilibria of F z,τ are strict saddles. Thus, following similar arguments as the proof of Theorem 3 in [42], we conclude that the set of all initial conditions that converge to these saddles has zero Lebesgue measure, which implies that almost every trajectory of the system converges to q d , i.e., the only stable equilibrium of (17), thus completing the proof.
Appendix A.7. Proof of Theorem 2 We begin by noting that, by virtue of (21), we only need to study the trajectories of (19) in the workspace's image, since that motion is traced exactly by our robot. Considering the first part of the Theorem 2, we note that by following the same arguments as in the proof of Proposition 3, we may conclude that the robot cannot escape throught the workspace's outer boundary. Likewise, assuming that q → q i for some i ∈ I N implies that (n θ T J T ∇ q ψ) T nθ T (q − q i ) cannot be identically zero in a neighborhood of q i . As such, sincek i ≥ 0 in the neighborhood of q i , k i cannot vanish as the robot approaches q i , which contradicts our original supposition.
To prove the second part of the Theorem 2, first we show that the only equilibria of the closed-loop system coincide with the critical points of ψ. Assuming that s v = 0, it is readily seen that both inner products in (20) vanish simultaneously only when ∇ q ψ = 0. Considering now the case when s v = 0, we note that this can only happen when q ∈ ∂D and nθ is tangent to ∂D. Forω to also vanish when s v = 0, the gradient ∇ q ψ should also be tangent to ∂D. Recalling that the adaptive laws for k ensure that ∇ q ψ will eventually point inwards, we conclude that no equilibria other than the critical points of ψ exist.
Next, we consider ψ as a lyapunov candidate function, whose derivative along the systems trajectories is given by (A6) (note that ψ does not depend on θ). Substituting (20) into the first term of (A6) yields: Regarding the remaining terms of (A6), given that g v,i ≤ nθ T ∇ q ψ 2 , one can readily verify that: Thus, invoking Lyapunov's Stability Theorem (Theorem 3.1 [41]) and LaSalle's Theorem (Theorem 3.4 [41]) concludes the proof similarly to Proposition 4.